rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata      = paste0(WD,"/_individual_data/_spfrawdata/us_spf_pgdp_ind.Rda","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(dplyr)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)
library(matrixStats)
theme_set(
  theme_minimal() +
    theme(legend.position = "right", axis.line = element_line(colour = "black"),text=element_text(size=12), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), )
)

ntile_na <- function(x,n)
{
  notna <- !is.na(x)
  out <- rep(NA_real_,length(x))
  out[notna] <- ntile(x[notna],n)
  return(out)
}


##### main text output  #######
load(indata)
nn              = 10

data$ID         = data$id
data            = data[data$year>=1970,]
data            = data[data$year<=2020,]

data = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))

min_id          = min(data$id)
max_id          = max(data$id)
fc_error_mean   = rep(NA,max_id-min_id+1)
fc_rev_mean     = rep(NA,max_id-min_id+1)
cons_mean       = rep(NA,max_id-min_id+1)
iter = 0
for (ii in min_id:max_id){
  iter   = iter + 1
  tmp    = data[data$id == ii,]
  fc_error_mean[iter] = mean(tmp$fc_err,na.rm=TRUE)
  fc_rev_mean[iter]   = mean(tmp$fc_rev,na.rm=TRUE)
  cons_mean[iter]     = mean(tmp$cons,na.rm=TRUE)
}
id              = min_id:max_id
df_tmp          = data.frame(id, fc_error_mean, fc_rev_mean, cons_mean)

data            = merge(data,df_tmp)

data$y       = data$fc_err - data$fc_error_mean 
data$x       = data$fc_rev - data$fc_rev_mean
data$z       = data$cons - data$cons_mean 

xbin = c(-20,-2,-1.5,-1.0,-0.5,0,0.5,1.0,1.5,20)
ybin = binMeans(data$y, data$x, bx =xbin, na.rm=TRUE)
xbin = c(-2.5,-1.75,-1.25,-0.75,-0.25,0.25,0.75,1.25,2.5)
df_tmp = data.frame(xbin, ybin)

pdf(file="_figures/figure_1_lhs.pdf")
plot(df_tmp$xbin, df_tmp$ybin, pch = 19, cex = 1.5, col = "black", main = "", xlab = "Revision", ylab = "", bty="l", cex.lab=1.2, xlim = c(-2.5,2.5), ylim = c(-0.4,0.4))
abline(lm(ybin ~ xbin, data = df_tmp), col="gray",  lwd=4)
dev.off()

#pdf(file="_figures/figure_1_lhs.pdf")
#plot(df_tmp$xbin, df_tmp$ybin, pch = 19, cex = 1.5, col = "darkblue", main = "", xlab = "Revision", ylab = "", bty="l", cex.lab=1.2, xlim = c(-2.5,2.5), ylim = c(-0.4,0.4))
#abline(lm(ybin ~ xbin, data = df_tmp), col="darkred",  lwd=4)
#dev.off()

zbin = c(-10,-2,-1.0,-0.5,0,0.5,1,2,3,4, 10)  
ybin2= binMeans(data$y, data$z, bx =zbin, na.rm=TRUE)
zbin = c(-1.75,-1.50,-0.75,-0.25,0.25,0.75,1.5,2.5,3.5,5.0)
df_tmp = data.frame(ybin2, zbin)


pdf(file="_figures/figure_1_rhs.pdf")
plot(df_tmp$zbin, df_tmp$ybin2, pch = 19, cex = 1.5, col = "black", main = "", xlab = "Consensus", ylab = "", bty="l", cex.lab=1.2, xlim = c(-2.5,5), ylim = c(-1.25,1.0))
abline(lm(ybin2 ~ zbin, data = df_tmp), col="gray",  lwd=4)
dev.off()

#pdf(file="_figures/figure_1_rhs.pdf")
#plot(df_tmp$zbin, df_tmp$ybin2, pch = 19, cex = 1.5, col = "darkblue", main = "", xlab = "Consensus", ylab = "", bty="l", cex.lab=1.2, xlim = c(-2.5,5), ylim = c(-1.25,1.0))
#abline(lm(ybin2 ~ zbin, data = df_tmp), col="darkred",  lwd=4)
#dev.off()
